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Abstract 

The Helmholtz equation arises in the study of electromagnetic radiation, optics, 
acoustics, etc. In spherical coordinates, its general solution can be written as a spherical 
harmonic series which satisfies the radiation condition at infinity, ensuring that the 
wave is outgoing. The boundary condition at infinity is hard to enforce with a finite 
element method since a suitable approximation needs to be made within reasonable 
distance from scatterers. Luckily, the Helmholtz equation can be represented as a 
Lippmann-Schwinger integral equation which removes the necessity of the boundary 
approximations and its Green's function can be expanded as a spherical harmonic 
series which leads to our numerical scheme based on spherical harmonic polynomial 
transform. In this paper, we present an efficient solver for the Helmholtz equation which 
costs 0(N log N) operations, where N is the number of the discretization points. We 
use the fast spherical harmonic transforms which are originally developed in [33J. The 
convergence order of the method is tied to the global regularity of the solution. At 
the lower end, it is second order accurate for discontinuous material properties. The 
order increases with increasing regularity leading to spectral convergence for globally 
smooth solutions. 

Keywords: Helmholtz equation, fast spherical harmonic transform, addition theorem, 
Lippmann-Schwinger integral equation, radiation condition, wave equation, 
Spectral convergence. 



1 Introduction 

Computational electromagnetics and acoustics are fundamental to understanding of many 
practical systems like scattering, microwave circuits, radar, antennas, remote sensing, seismic 
exploration, ultrasound and tomography. Therefore, the demand for efficient numerical sim- 
ulations of electromagnetic and acoustic fields is only increasing. The electromagnetic field is 
the solution of Maxwell's equations which are coupled with more than one unknown. These 
equations can be uncoupled by raising their order resulting in the wave equation. Thus, 
the Helmholtz equation reduced from the wave equation has been a classical problem ([5], 
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[23] ) to solve and central focus of research for many decades and still drawing ongoing atten- 
tions. Using the separation of variables, the Helmholtz equation has been analytically solved 
for a simple geometry with homogeneous medium. In particular, for a spherical geometry, 
the solution is analytically given as the spherical harmonic series expansion in the angular 
direction and Bessel series expansion along the radial direction ([5]) which inspires our nu- 
merical method based on spherical harmonic series expansion that is particularly efficient for 
spherical objects. There have been many efforts to develop fast polynomial transformations 
([12], [15], [22], [19], [33]) to have equivalent speed advantages like FFTs. Here, we used the 
fast spherical harmonic transforms ([19]) developed in [33] which is based on fast multipole 
method ([15]) to have 0(N log N) costs comparable to those based on FFTs while enjoying 
that the separation of variables that results from the addition theorem ([ID]). This readily 
translates into quadratures that converge with higher order than those implicit in FFT-based 
schemes. 

For the Helmholtz equation concerned in this paper is to find the scattered field generated 
due to the incident field with outgoing radiation condition. The most popular numerical 
methods are finite element methods (FEM) ([13], [18], [21], [26], [31]) and integral equation 
methods (IEM) ([I], [7J, [8], [19], [20], [2S], [35]). Although FEM can handle arbitrarily 
shaped obstacles with ease and the resulting matrix is sparse, it imposes a serious difficulty 
in solving scattering problems arising from the infinite size of domain. For this reason, 
great effort has gone into the design of approximate local boundary conditions and perfectly 
matched layers ([2], [3], p2]) that minimize spurious reflections which is by no means a 
trivial matter Q18J). Integral equation methods (IEM), on the other hand, implicitly ac- 
count for radiation conditions through the use of outgoing Green's functions. This very use 
of singular Green's functions, on the other hand, also translates into numerical challenges. 
Moreover, these methods lead to a linear system involving a full matrix and thus they are 
not competitive unless a specialized strategy is used to accelerate matrix-vector products; 
examples of accelerated IEM include those based on FFTs ([7], [20], [28J) and those that use 
fast multipole expansions ([IB]. [25] . [27]). 

In this paper, inspired by the work in ([6], [7J, |20j), we present a new accelerated IEM based 
on the addition theorem (PU]) and fast spherical harmonic transforms ([IS], [S3])- The con- 
vergence rate of our new algorithms is tied to the global regularity of fields. In particular, 
they converge with second order for the most singular case of discontinuous material prop- 
erties and with increased rates for more regular arrangements; for smooth configurations 
the convergence is spectral. This paper is organized as follows. The Lippmann-Schwinger 
integral equation is presented in Section 2. The numerical factorization of the integral equa- 
tion based on spherical harmonic series expansion is explained in Section 3. In Section 4, 
numerical implementations and their expected costs are derived, and in Section 5, numeri- 
cal examples are given to confirm the predicted performance of the algorithms described in 
Section 4. Finally, in Section 6, the content of this paper is summarized. 
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2 Lippmann-Schwinger integral equation 

The wave equation states that 

where n(r) is the refractive index and Co is the propagation speed of the wave in air. If we 
assume that the wave function is time-harmonic as 

u(r,t) = u(r)e- iwt , 

then a spatial solution u(r) satisfies the Helmholtz equation 

Au + k 2 n 2 u = in M 3 . (2) 

Here we will consider the scattering problem to determine the total field generated by a 
given incident field u % . Then the total field u which is the sum of u % and the scattered field 
u s satisfies equation (j2J) while u l is the solution of 

Am* + fcV = in M 3 . (3) 

The Sommerfeld radiation condition is given at infinity which guarantees that the scattered 
field is outgoing, 

f du s \ 

lim r — iku s = 0. (4) 

r^oo \ or J 

This boundary condition must be approximated at the reasonable distance from scatterers 
which is the main challenge to solve the Helmholtz equation. But, this problem can be 
avoided if one appeals to the equivalent Lippmann-Schwinger integral equation which states 

u(x) = u\x) - k 2 J y)u(y)m(y) dy, x E M 3 (5) 
h 

where m = 1 — n(x) 2 and the inhomogeneity is of compact support Q. Therefore, we will 
develop an efficient solver based on the Lippmann-Schwinger integral equation. 



3 Spherical harmonic series expansion 

To solve the integral equation ((Sj), we resort to the spherical harmonic series expansions of 
the Green's function The addition theorem in pU] states that 

1 e ik\x-y\ 

$(x,y) = — 1 

47r | x — y | 

oo n (6) 
n=0 m=—n 
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where 

x = p(sin 9 cos <p, sin 9 sin <p, cos 6) , 
y = p(sin 6 cos <p, sin 9 sin <p, cos 9) , 
p< = min(p, p) and p> = max(p, p). 
So, if we approximate w and u % by a truncated spherical harmonic series as 

F n 

u F ( P , M} = EE <iP)Yn(0, <P) (7) 

n=0 m=— n 
F n 



u^( P ,9,<p)} = j2 E < m (p)^ m 



71 V." ) 

n=0 m=— n 



then from equation ([5]) and the orthogonality properties of the spherical harmonics, m{r) 
can be approximated without losing accuracy as 



2F ri 

,2F 



( P ,M)} = E E <(p)WM)- ( 9 ) 



m 

n=0 m=— n 



Therefore, the formulation in ([5]) becomes 

<(p)=< m (p)+^n m (p), (10) 



where 



R 



K?(p) = -k 3 J h£\k P> )j n (k P< )I™( P ) P 2 dp, (11) 



and 



2-7T 7T 



C(p)= y y ^(p,^,0)m^(p,^,0)y n m (^,0)sin^^0, (12) 

0=0 0=0 

n = 0, 1, . . . , F. To solve the formulation in ffTQ|) . we use the linear solver GMRES which 
requires the fast evaluation of angular integration of I™(p) in f|T2|) and the radial evaluation 
of K™(p) in <HH). 

4 Numerical implementation 

4.1 Angular integration 

Due to the orthogonality of spherical harmonics, equation ([12]) can be written as 

3F n 

u F (p, 9, 0)m 2F (p, 9, 0) = E E COO W (13) 



Therefore, we define spherical harmonic transform as 



{m,<PM=0 ^ {{Cn} n m =-n}Lo (14) 



'31 i 1,3 

where 



F n 



n=0 m=—n 



and its inverse, for appropriate choices of the angles {6^,0^}; see [33] , [T9] . 
In [33], it is shown that if we define ST as 



then 

AT-1 

^C^ +n (xf), je {0,l,2,,iV-l} (15) 



n=0 

can be computed in 0(N log N) operations for arbitrary xf. 

This work leads to fast spherical harmonic transform (FSHT) and its inverse (IFSHT) and 
these polynomial transforms cost 0(F 2 logF) for (TTJJ. Thus, /™(p) of f[T2|) is computed as 
follows 



{{C(p)}^=-Jn=o = FSHT 3F (IFSHT 3F ({{u™(p)r m= _X=o 
■ IFSHT 3F ({{m^(p)} n m= _ n } 2 n F =0 )). 



( 16 ) 



Therefore, the angular integration costs 0(F 2 logF) for each p. 
4.2 Radial integration 

The radial integral K™(p) in (TTTT) has a corner-type singularity at p = p therefore, we write 
it as 

min(a,R) R, 

Z ^P^ = h^\ka) J Jn (kp)C(p)p 2 dp + Jn (ka) J hH\kp)C(p)p 2 dp 

min(a,R) 
min(a,R) R 

min(a,R,) 



y n (fca) y j n {kp)I™(p)p 2 dp + j n (ka) J y n {kp)I™(p)p 2 dp 



R 

m/ \ 2 



+ Jn(fca) Jn (kp)I™(p)p 2 dp. 







Although it is numerically challenging to evaluate the Hankel function hn\kp), the product 
of the jn(kp) and hn\kp) is bounded which leads us to define modified Bessel functions 



5 



j n (p) and y n (p) in the form 
^ l-3-5-...(2n + l) , 

JAP) ■= MP) 



p fl+1 



¥ 2 



+ 



(|P 2 ) 2 



l!(2n + 3) 2!(2n + 3)(2n + 5) 

(|P 2 ) 2 



+ . . . 



1 - 



+ 



l!(l-2n) 2!(l-2n)(3-2ra) 



+ 



and with these modified Bessel functions we obtain 

min(a,R) 



K?(a) = 

nK J 2n + l 



y~(ka) J {^j n+l j n {kp)ln(p)k 2 pdp 



(18) 



R 



+ jn(ka) 



min(a,R) 



(-) n y n (kp)I™(p)k 2 pdp 



R 



+ j n (ka)(ka) 



{kp) n {-2n - l)k 3 
l-3 2 -5 2 ...(2n+l) 3 



j n (kp)i™( P yd P . 



(19) 



For the radial integration, we divide the integration domain in a number TV, of equi-length 
interpolation intervals Uj = uj], 1 < j < N t on which we approximate for p E Uj, 

N d -1 

where 



(20) 



Tf" U '(A>)=Ii 



p-(^+^ )/2 
' l («}-«J)/2 



are the Chebyshev polynomials in Uj. It costs 0(N d (logN d )N i ) to compute the Chebyshev 
coefficients of I™(p) for fixed n and m. If we denote {p^}^ the Chebyshev points in Uj, 
the following equation holds 



j p n+2 Jn(kp)C(p)d P = J p n+2 Ukp)I™{p)d P 





N d -1 



+ 



Y.4,m,n I P n+2 Ukp)T^ U \p)dp. 

1 — n ** 



(21) 



1=0 



Therefore, if the moments J p n+2 j n (kp)T l 3 ' 3 (p) are computed and stored, then it costs 



pi 



pi 



O(NjNi) to compute / p n+2 j n (kp)I™(p) dp for 1 < j < N t and 1 < k < N d . Similarly, the 
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R 

computation of j p 1 ~ n y n (kp)I™(p) dp also costs O(NjNi) for 1 < j < iVj and 1 < k < N d . 

Pi 

With these moments and precomputed coefficients Ci(a), C 2 (a) and C 3 (a), K™(a) can be 
written as 

min(a,R) 

K™(a) = C 1 (a) j p n+2 Ukp)In{p)dp 
o 

R 

1-n 



+ C 2 (a) / p 1 - n y n (kp)I^(p)dp (22) 



min(a,R) 
R 

„n+2 



+ C 3 (a) j p n+2 j n (kp)I™(p)dp. 

Q 

Therefore the total cost for radial integration is 

0[iV i JV d F 2 (logiV ci + iV < , + l)] =0(JV). 

5 Numerical examples 

To show the predicted performance of the algorithm, four examples are presented below. 

Example 5.1 Consider the scattering off a homogeneous sphere of radius 1 and the refrac- 
tive index n(x) = 2. In this case, the problem is explicitly solvable, therefore comparison 
with the exact solution is possible. A plane wave incident in the positive ^-direction can be 
written as 

oo 

e ifc*-(o,o,i) = e ik P cos8 = J2i n (2n + l) Jn (kp)P n (cos6). (23) 

n=0 

If we differentiate equation ( )23|) rrii nc times with respect to t = cos 8, we obtain the incident 
wave extended from ([23]) as 



yi _ p\m inc \ e ikp cos 6 s ^\m inc \ Q^,m inc 4> 

_ i"(2n + 1) V4tt(n + |m„, c |)! (24) 

For this incidence the exact solution [5] is 

'Er=im inc i{^n(2fcp) -q^mwmm), p < 1 



TOO 



where 



Figure 1: The incident field intensity \u l \ 2 where m inc = 1 for examples 5.1-5.3. 



By enforcing C 1 continuity of u s along the material discontinuity, {a n } and {&„} are deter- 
mined. The numerical error is computed as E = \\u exact — M ap prox ||oo between the approximate 
and exact solutions for m inc = 1 and k = 5 in fl24"l) and different values of the interpolation 
orders N<i in (1201) in Tables 1-3. 



Ni 


time per iteration 


GMRES iteration 


relative error 


error ratio 


2 3 


54 (sec) 


34 


0.564482 




2 4 


108 (sec) 


35 


0.20477 


2.81207 


2 5 


217 (sec) 


36 


0.0532706 


3.7178 



Table 1: Radial convergence for Example 5.1: the sphere centered at the origin. 
Parameters: m inc = 1, k = 5, F = 2 5 - 1, N d = 2, < p < 2, 
GMRES tolerance = le - 10. 





time per iteration 


GMRES iteration 


relative error 


error ratio 


2 3 


108 (sec) 


35 


0.00192818 




2 4 


217 (sec) 


36 


0.000193606 


10.0197 


2 5 


435 (sec) 


36 


1.35701e-05 


14.1395 



Table 2: Radial convergence for Example 5.1 for = 4, GMRES tolerance = le — 10, 
Same parameters as in table 1. 



Figure 2: The field intensity \u\ 2 for example 5.1 where m inc = 1, k = 5. 





time per iteration 


GMRES iteration 


relative error 


error ratio 


2 3 


218 (sec) 


69 


2.92336e-08 




2 4 


435 (sec) 


67 


1.70495e-10 


171.523 


2 5 


1098 (sec) 


66 


6.61144e-13 


257.404 



Table 3: Radial convergence for Example 5.1 for = 8, GMRES tolerance = le — 15, 
Same parameters as in table 1. 



Example 5.2 To test the convergence when there is a material discontinuity for the angular 
integration, we consider the same sphere as in example 5.1 but now centered at (0, 0, d). 
For this case, we can obtain the incident wave as 



and the exact spherical harmonic series expansion for m(x) as follows 
m(p,9) 



n=l v ' 

1 



+ -(l-cos# )v / 47ry °(fl, 



(25) 



where < (d — r) < p < (d + r), n Q = 2 is the refractive index and (psin# , pcos^o) is a 
solution of 

2,2 2 2 i / j\2 2 

y +z =p,y +{z-d) =r . 
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The exact solution is obtained from the spherical harmonic transform of the exact values 
computed from the solution in example 5.1 with shifting. In tables 4-6 we present radial 
and angular convergence studies for this case where the discontinuity does not lie on the grid. 





time per iteration 


GMRES iteration 


relative error 


error ratio 


2 5 


1203 (sec) 


26 


0.299451 




2 6 


2408 (sec) 


27 


0.0817522 


3.64575 


2 7 


4987 (sec) 


27 


0.020442 


3.99036 



Table 4: Radial convergence for Example 5.2: the sphere centered at (0,0,2) with radius 1. 
Parameters: m inc = 1, k = 5, F = 2 7 - 1, N d = 2, < p < 4, 
GMRES tolerance = le - 5. 





time per iteration 


GMRES iteration 


relative error 


error ratio 


2 4 


601 (sec) 


6 


2.88611e-05 




2 5 


1203 (sec) 


6 


8.11132e-06 


3.51951 


2 6 


2406 (sec) 


6 


1.97643e-06 


4.0819 



Table 5: Radial convergence for Example 5.2: the sphere centered at (0,0,2) with radius 1. 
Parameters: m inc = 3, k = 1, F = 2 7 - 1, N d = 2, < p < 4, 
GMRES tolerance = le - 10. 



F 


time per iteration 


GMRES iteration 


relative error 


error ratio 


2 4 -l 


352 (sec) 


21 


1.9425 




2 5 -l 


866 (sec) 


26 


0.113651 


17.0918 


2 6 -l 


2061 (sec) 


27 


0.00157294 


72.2536 



Table 6: Angular convergence for Example 5.2: the sphere centered at (0,0,2) with radius 1. 
Parameters: rrii nc = 1, k = 5, N d = 2, iVj = 2 7 , < p < 4, 
GMRES tolerance = le - 5. 



Example 5.3 To test the convergence of a non-spherical object, we consider a square rotated 
by 45 degree and axisymmetric along z direction. The refractive index n(r) is 2 in + < 1 
shown in Figure 3. The incident wave is the same as in example 5.1. We present radial and 
angular convergence studies in tables 7-9. 
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Figure 3: The scatterer for example 5.3. 




Figure 4: The field intensity \u\ 2 for example 5.3 where irii nc = 1, k = 5. 
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time per iteration 


GMRES iteration 


relative error 


error ratio 


2 4 


259 (sec) 


20 


0.606676 




2 5 


559 (sec) 


20 


0.0568001 


10.6809 


2 6 


1081 (sec) 


20 


0.0416413 


1.36403 


2 7 


2072 (sec) 


20 


0.0126427 


3.29371 



Table 7: Radial convergence for Example 5.3: the square rotated by 45 degree. 
Parameters: m inc = 1, k = 5, F = 2 6 - 1, N d = 2, < p < 2, 
GMRES tolerance = le - 5. 





time per iteration 


GMRES iteration 


relative error 


error ratio 


2 7 


2072 (sec) 


8 


0.000389159 




2 8 


4135 (sec) 


8 


0.000123866 


3.14176 


2 9 


8729 (sec) 


8 


1.79324e-05 


6.90739 


2 io 


17823 (sec) 


8 


3.94232e-06 


4.54869 



Table 8: Radial convergence for Example 5.3: the square rotated by 45 degree. 
Parameters: m inc = 1, k = 1, F = 2 6 - 1, N d = 2, < p < 2, 
GMRES tolerance = le - 10. 



F 


time per iteration 


GMRES iteration 


relative error 


error ratio 


2 4 -l 


177 (sec) 


26 


0.0328253 




2 5 -l 


436 (sec) 


26 


0.00439355 


7.47126 


2 6 -l 


1036 (sec) 


26 


0.000548512 


8.00994 



Table 9: Angular convergence for Example 5.3: the square rotated by 45 degree. 
Parameters: m inc = 1, k = 5, N d = 2, N { = 2 6 , < p < 2, 
GMRES tolerance = le - 10. 



Example 5.4 The final example is to show the dependence of the scheme on the regularity 
of the scattering medium. We consider a refractive index given by 

. f (1+ I costf f sin 1 ™-/ 1 # e im "^) 1/2 , 1 < p < 2 
n(p,cos9) = < (26) 
II, < p < 2 or p > 2. 



Then, 

-|f|^(l-* 2 ) i!= T Zi e <r, w*, l<p<2 



. . ,. I m refl 

! 1 / ; >i \ _ t-- \ 

m(p,t) 



0, 0<p<lorp>2 
12 



and the exact spherical harmonic series expansion for m(p, t) is given by Yl'iLo m2 i(p)Y\^ ref f \ + 2i(t)i 
where 



/2-7T /»7T 
/ Y 2 ref '\+2i( e ' <t>) cos/3 9 sin |m '' e/l 0e iro ™'* cos 9 dOdcf) 
r =oJe=o re 



(2n + l)(n - 


m re/ 


)! 7T + 1 


(n + 


m re/ |)! 2"+l"wl r(l + f) r(| + f) 



(27) 



(n + |m re/ |)! 



ni=o /l+ ' _1 (^+i+«) ( W "Ke/D ! ' 



In tables 10-12, we present the order of convergence for (3 =0.4, 1.4, 2.4 to show the corre- 
lation between smoothness of the refractive index n(x) and the order of convergence. 



F 


time per iteration 


GMRES iteration 


relative error 


log 2 (error ratio) 


2 3 -l 


83 (sec) 


6 


0.00186872 




2 4 -l 


419 (sec) 


6 


7.42083e-06 


7.97625 


2 5 -l 


2048 (sec) 


6 


1.19402e-06 


2.63575 



Table 10: Angular convergence for Example 5.4 for (3 = 0.4 (m G C 0A , u G C 2A ). 

Parameters: m inc = 3, m re f = 1, k = 0.5, Na = 8, Ni = 4, < p < 4, 
GMRES tolerance = le - 10. 



F 


GMRES iteration 


relative error 


log 2 (error ratio) 


2 3 -l 


5 


0.0018687 




2 4 -l 


5 


6.99546e-07 


11.3833 


2 5 -l 


5 


5.82286e-08 


3.58662 



Table 11: Angular convergence for Example 5.4 for (5 = 1.4 (m G C 1A , u G C 34 ). 
Same parameters as in table 10. 



F 


GMRES iteration 


relative error 


log 2 (error ratio) 


2 3 -l 


5 


0.00186869 




2 4 -l 


5 


6.6957e-08 


14.7684 


2 5 -l 


5 


2.74242e-09 


4.60971 



Table 12: Angular convergence for Example 5.4 for (3 = 2.4 (m G C 2A , u G C AA ). 
Same parameters as in table 10. 
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Figure 5: The field intensity \u\ 2 for example 5.4 where rrii nc = 1, m re f = 7, — 0.2, k — B. 

6 Summary 

In this paper, an efficient solver for scattering by penetrable three-dimensional structures 
is presented. The solution is obtained by the iterative evaluation of Lippmann-Schwinger 
integral equation and its efficiency comes from the use of the addition theorem and fast 
spherical harmonics transforms. The scheme allows for such evaluations in 0(N log N) op- 
erations, where N is the number of the discretization points. The convergence order of the 
method, on the other hand, is tied to the global regularity of the solution. At the lower 
end, it is second order accurate for discontinuous material properties. The order increases 
with increasing regularity of the refractive index leading to spectral convergence for globally 
smooth solutions. 
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